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ORIGIN OF THE COSMIC RAY SPECTRAL HARDENING 
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ABSTRACT 

Recent data from ATIC, CREAM and PAMELA indicate that the cosmic ray energy spectra of protons and 
nuclei exhibit a remarkable hardening at energies above 100 GeV per nucleon. We propose that the hardening 
is an interstellar propagation effect that originates from a spatial change of the cosmic -ray transport properties 
in different regions of the Galaxy. The key hypothesis is that the diffusion coefficient is not separable into 
energy and space variables as usually assumed. Under this scenario, we can reproduce the observational data 
well. Our model has several implications for the cosmic-ray acceleration/propagation physics and can be tested 
by ongoing experiments such as the Alpha Magnetic Spectrometer or Fermi/LXT, 
Subject headings: acceleration of particles — cosmic rays — diffusion — turbulence 



1. INTRODUCTION 

Understanding the spectral features in cosmic rays (CRs) is 
of fundamental importance for any theory of their origin and 
propagation. The spectrum of CR nuclei above ~ 10 GeV per 
nucleon is thought to be the result of acceleration mechanisms 
in supernova remnants (SNR), which is steepened by propa- 
gation in the interstellar medium (ISM) due to leakage from 
the Galaxy (Strong et al. 2007). The acceleration of primary 
CRs is described by the diffusive shock acceleration (DSA) 
theory that predicts power-law spectra Q oc E~ v with slope 
v k, 2.0-2.2. The CR transport in the turbulent magnetic fields 
is modeled as a diffusion process, with coefficient K oc E s , 
into a magnetic halo of typical height 2L ~ lOkpc. Inter- 
actions of CRs with the ISM give rise to secondary nuclei 
such as Li-Be-B. This description predicts power-law spec- 
tra such as ~E~ V ~ S for primary nuclei {e.g., H, He) and 
^E~ s for secondary-to-primary ratios (e.g., B/C). The pa- 
rameters v and S are not firmly predicted a priori and can be 
inferred from the data. Observations constrain v+5 ps 2.7 (de- 
pending on the element) and 5 — 0.3-0.7, whereby v — 2.0- 
2.4. Yet, this picture dramatically over-predicts the CR 
anisotropy at >TeV energies, which would suggest an al- 
most energy-independent diffusivity. Although the CR flux 
may also be affected by stochasticity effects on SNR events 
(Blasi & Amato 2012b), the high degree of isotropy observed 
in CRs poses a serious challenge to the conventional approach 
to the CR propagation. Even more dramatically, the CR spec- 
tra at E > lOOGeV/nucleon exhibit a remarkable hardening 
with increasing energies. This feature has been established 
by recent experiments ATIC-2 (Panov et al. 2009), CREAM 
(Ahn et al. 2010) and PAMELA (Adriani et al. 201 1), though 
the data are a little discrepant with each other. The proposed 
explanations interpret the hardening as a source effect con- 
nected with acceleration mechanisms (Biermann et al. 2010), 
nearby SNRs (Thoudam & Horandel 2012) or arising by dif- 
ferent populations of CR sources (Zatsepin & Sokolskaya 
2006; Yuan et al. 201 1). The subsequent impact of the hard- 
ening on the secondary CR production was estimated by 
Lavalle (201 1) for e+, Donato & Serpico (201 1) for p and 7- 
rays, and Vladimirov et al. (201 1) for light nuclei. The latter 
also considered a possible diffusive origin of the effect, de- 
scribed by means of an effective break in the slope of K(E), 



around -200 GeV/nucleon, from 5 = 0.3 to 5 = 0.15. 

Contrary to these works, we propose that the hardening 
originates from a spatial change of the CR diffusion prop- 
erties in the different regions of the Galaxy. In fact, propa- 
gation studies assume that CRs experience the same type of 
diffusion in the whole propagation region, with only few ex- 
ceptions for particular environments (e.g., in SNR shells or 
inside the heliosphere) that are treated separately. A more 
realistic position-dependent diffusivity (correlated with the 
SNR density) is considered in some works (e.g., Shibata et al. 
2004; Gebauer & de Boer 2009; Di Bernardo et al. 2010), but 
they still adopt a unique energy dependence for K in the 
whole halo, which leads to essentially unchanged results for 
the shapes of CR spectra at Earth. Such descriptions may 
represent an over-simplification of the problem. In theoret- 
ical considerations, the diffusion is caused by the CR scat- 
tering on hydromagnetic waves which, in turn, depends on 
the nature and scale distribution of the magnetic -field irreg- 
ularities. While SNR explosions may generate large irregu- 
larities in the region near the Galactic plane, the situation in 
the outer halo is much different because there are no SNRs. 
The main source of turbulent motion in the halo is presum- 
ably represented by CRs themselves. From these considera- 
tions, Erlykin & Wolfendale (2002) found that the turbulence 
spectrum in the halo (in the far outer Galaxy) should be flat- 
ter than that in the plane (in the inner Galaxy). This implies 
strong latitudinal changes (gradual radial changes) for the pa- 
rameter S and suggests spectral variations of CRs throughout 
the Galaxy. Noticeably, new data reported by Fermi/LAT on 
the diffuse 7-ray emission at — 10-100 GeV of energy seem 
to support these suggestions. The 7-ray spectra observed 
near the Galactic plane (latitude |6| <8°) are found to be 
harder than those at higher latitudes, and similar differences 
are also found between the inner Galaxy (longitude I < 80° 
or I > 280°) and the outer Galaxy (Ackermann et al. 2012). 

In this Letter, we focus our attention on the latitudinal 
changes of the CR diffusion properties (which are expected 
to be more extreme) and we show that they inevitably lead 
to a pronounced hardening for the energy spectra of CR nu- 
clei at Earth. Using analytical calculations, we examine the 
implications of this scenario for the main CR observables - 
primary spectra, secondary-to-primary ratios and latitudinal 
anisotropy- and discuss its connections with the open prob- 
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lems in the CR acceleration/propagation physics. 

2. CALCULATIONS 

We use a simple model of CR diffusion and nuclear inter- 
actions. The effects of energy changes and convection are 
disregarded. The Galaxy is modeled to be a disk, with half- 
thickness h, containing the interstellar gas (number density n) 
and the CR sources. The disk is surrounded by a diffusive halo 
of half-thickness L and zero matter density. For simplicity we 
give a one-dimensional description (infinite disk radius) in the 
thin disk limit (h <C L) and we restrict to stable species. For 
each CR nucleus, the transport equation reads 



dN 
~dt 



|- (k(z)^] -2hS(z)T mcl N + 2h6(z)Q, (1) 
oz \ oz I 



where N(z) is its number density as function of the z- 
coordinate, K(z) is the position-dependent diffusion coeffi- 
cient, and r mcl = f3cna lncl is the destruction rate in the ISM 
at velocity /3c and cross section a 1 . The source term Q 
can be split into a primary term Q pr i, from SNRs, and a sec- 
ondary production term Q scc = ^\ Tj' 3 ' Nj, from spallation 
of heavier (j) nuclei with rate r® pa . The quantities N, K, 

Q, and r incl depend on energy too. Since no energy changes 
are considered, such a dependence is only implicit and can be 
ignored for the moment. To solve Eq. 1 we assume steady- 
state conditions (dN/dt = 0). We define ai(z) = K'/K, 
a 2 = —2hT lnel /K , and a 3 = 2hQ/K , where we have de- 
noted K' = dK/dz and K Q = K(z=0). In the halo (z § 0) 
the equation reads N" + a.\ N' = 0, which is readily solved 
as N±(z) = p± + u±X(z), where the subscripts =p indicate 
the solutions in the z ^ half-planes. The function X(z) is 
defined as 



X(z) 



So a 1 (z" )dz" dz > 



(2) 



The boundary conditions, N(±L) = 0, provide the relation 
u± = —p/A±, where A± = A(±L). From the continuity 
condition in the disk, iV_(0) = N + (0), one obtains p + = 
P- = p. Assuming that K(z) is an even function, one can 

see that X(z) must be odd. Thus we define X(z) = X(\z\) and 
A = X(L), so that u± = ^fp/A. To determine p, we integrate 
the transport equation in a thin region across the disk 

N' + (e)-N'_(-e) + J otiN'dz + a 2 N(0) + a 3 = . (3) 

The limit e — > gives p = az/ (2A" 1 — 0.2). After replacing 
ot\, c*2, and a 3 with the original quantities, the solution reads 



N(z) 



Q 



En. 

hA 



1 



X(z) 



A 



The function X(z) can be expressed as 

=1 



X(z) 



"J? K '/K**dz' = K Q 



dz 
K{z) 



(4) 



(5) 



From this toy model, one can recover the homogeneous 
diffusion model (HM) by setting K' = 0, which gives A = \z\ 
and A = L. Simple models of inhomogeneous diffusion 
can be described by a diffusion coefficient of the type 
K(z,E) = J{z)Kq(E). For example, Di Bernardo et al. 



(2010) adopt f(z) = e) z \l Zt . In this case one finds 
A = Zt (l — e - ^'/ 2 *), where the limit z t 3> L recovers 
the HM, and the limit z t <C L gives A w z t . The latter 
case provides a more natural description of the latitudinal 
CR density profile, as it is insensitive to the halo bound- 
aries ±L. However, the model predictions in terms of CR 
spectra at Earth remain equivalent to those of the HM after 
a proper choice of A. This is a general property of Eq. 5: 
as long as K(z,E) is separable in z and E, the function 
A is independent on energy and the spectra at z=0 are 
equivalent to a mere rescaling of the model parameters. 
We remark that the energy-space variable separation is 
implicitly assumed in all CR propagation models. Physi- 
cally, it describes a unique diffusion regime in the whole 
halo, given by Kq(E), while f(z) allows for spatial varia- 
tions in its normalization. The quantity A can be regarded 
as an effective halo height experienced by CRs at equilibrium. 

We now put into practice our hypothesis on the latitudinal 
variations of the CR diffusion properties, which implies that 
K(z, E) is not separable as f{z)K$(E) everywhere. We fol- 
low the arguments given in Erlykin & Wolfendale (2002), but 
our aim is not to inspect the astrophysical plausibility of their 
suggestions. Rather, we consider a phenomenological sce- 
nario in order to illustrate the effect and its consequences for 
the main CR observables. We adopt a simple two-halo model 
(THM) consisting in two diffusive zones. The inner halo, rep- 
resenting the region influenced by SNRs, is taken to surround 
the disk for a typical size £L of a few hundred pc (£ ~ 0.1). 
The outer halo, representing a wider region where the turbu- 
lence is driven by CRs, is defined by £L < \z\ < L. The 
diffusion coefficient is taken of the type 



K(z,p) 



k (3 P s 



for \z\ < [inner halo) 
for \z\ > £L (outer halo) , 



(6) 



where p = R/Ro and ko specifies its normalization at the ref- 
erence rigidity Rq. We use the same fco values on both halos, 
as their relative normalization can be controlled by Rq. For 
R > Rq, the diffusion coefficient of Eq. 6 produces a higher 
CR confinement in the inner halo (with 6 ~ 1/3), whereas the 
outer halo (with A ^0.5 - 1) represents a reservoir from 
which CRs leak out rapidly and can re-enter the inner halo. 
From Eq. 5, one can compute A = X(z,p) and A = A(p). 
The latter reads 



A(p) = L[£+(l-0p 



-Al 



(7) 



That is, the effective halo height is a rigidity-dependent quan- 
tity that affects the model predictions at z = 0. This effect can 
be better understood if one neglects the term r mcl and takes a 
source term Q pl [ ~ p~ v . From Eq. 4 and Eq. 7 one finds 

No = N(z = Q)~± + (1 - H)p-"- S - A } , 

Kq 

(8) 

which describes the CR spectrum as a result of two compo- 
nents. Its differential log-slope as a function of rigidity reads 



= — n ~ v + d H - 

dlogp 1. 



A ' 



(9) 



which indicates a clear transition between two regimes. In 
practice, the low-energy regime (7 w v + 5 + A) is never 
reached due to spallation (neglected in the above equations) 
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that becomes relevant and even dominant over escape (r incl > 
In this case the log-slope is better approximated by 7 ~ 

v + i(5 + A) (Blasi & Amato 2012a). The hard high-energy 
regime (7 w v + 5) is determined by the diffusion properties 
of the inner halo only, because the outer halo is characterized 
by a much fast particle leakage. In this limit one has A ps £L. 
The effect vanishes at all rigidities when passing to the HM 
limit of £ — > 1 (one-halo) or A —> (identical halos), where 
one recovers the usual relation 7 = v + S. Furthermore, from 
Eq. 4, it can be seen that the intensity of the harder component 
diminishes gradually with increasing \z\, i.e., the CR spectra 
at high energies are steeper in the outer halo. 

3. RESULTS AND DISCUSSIONS 

We compute the THM spectra at Earth by 
J{E) = ^Nq(E), from Eq. 4, at kinetic energies 
above 10 GeV. The SNR energy spectra are taken as 
Qpri = Yp- 1 (R/R )- ,/ e- R / R ™°* , where R max rep- 
resents the maximum acceleration rigidity attainable by 
SNRs. The constants Y are determined from the data at 
~ lOOGeV/nucleon. The indices v are taken as Z-dependent 
to account for the observed discrepancies among elements. 
MalkovefaL (2012) and Ohira & Ioka (2011) gave strong 
arguments for ascribing such discrepancies to SNRs. The 
ISM surface density is taken as h x n = lOOpc x 1 cm -3 . 
The two halos are defined by L = 5kpc and £L = 0.5kpc, 
but the physical parameters that enter the model are ko/L 
and £, where both quantities are also degenerated with Rq. 
Concerning the diffusion parameters, we consider two case 
studies: THMa (5 =1/3 and A =0.55) which adopts a 
Kolmogorov-type diffusion in the inner halo, and THMb 



TABLE 1 
Model parameters 



parameters 


THMa 


THMb 


iv(H; He; CNO; Fe) 


2.29; 2.17; 2.17; 2.20 


2.43; 2.31; 2.31; 2.34 


Ro', i?max 


2GV; 2.510 6 GV 


2GV; 2.510 6 GV 


k /L 


0.007 kpc/Myr 


0.010 kpc/Myr 


8; A; £ 


1/3; 0.55; 0.1 


1/6; 0.55; 0.1 



(S = 1/6 and A = 0.55) to test the extreme case of a very 
slow diffusivity. The main parameters are summarized in 
Tab. 1. The practical implementation of the model follows 
Tomassetti & Donato (2012), see also Maurin et al. (2001). 

The energy spectra of H, He, CNO and Fe are shown in 
Fig. 1 in comparison with the data. Results are shown for 
THMb only (THMa and THMb predictions are indistinguish- 
able for primary CRs). Our calculations (solid lines) are in 
good agreement with the data within their uncertainties. At 
low energies (< lOGeV/nucleon) the solar modulation is ap- 
parent and it is described using a force-field modulation po- 
tential 4> =400MV (Gleeson & Axford 1968). Note, how- 
ever, that our model may be inadequate in this energy region 
due to the low-energy approximations. At higher energies, 
our model reproduces well the observed changes in slope, in 
agreement with the trends indicated by the data. It should be 
noted, however, that the sharp spectral structures suggested 
by the PAMELA data at ~ 300 GeV cannot be recovered. The 
THM predictions are also compared with HM power-law ex- 
trapolations (dashed line) to better illustrate the differences. It 
can be seen that the spectral upturn is slightly less pronounced 
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FIG. 2. — B/C ratio (a) and latitudinal anisotropy (b) for THMa, THMb 
and HM. B/C data are from HEA03-C2 (Engelmann et al. 1990), CREAM 
(Ahn et al. 2008), AMS-01 (Aguilar et al. 2010), ATIC-2 (Panov et al. 2007), 
and TRACER (Obermeier et al. 201 1). Air Shower and Underground fj, data 
are from Amenomori et al. (2005); SMP data are from Strong et al. (2007). 



for elements with large mass M such as Fe, due to the com- 
peting action of the term p ncl below a few lOOGeV/nucleon 
(note that T incl oc M 3 ). Above - lOOTeV/nucleon, the el- 
emental knees are matched well using i? max = 2.5T0 6 GV, 
which is somewhat lower than that from HM-based estimates 
(Blasi & Amato 2012a). For instance, the HM spectra of 
Fig. 1 (dashed lines) employ i? max =7-10 6 GV. 

A distinctive feature of our model is provided by the 
secondary-to-primary ratios. From Eq. 4, neglecting P , it 
is easy to see that, roughly, secondary-to-primary ratios must 
harden as oc A/Kq. In Fig. 2a we plot the B/C ratio. Re- 
markably, THMb describes well the data down to energies 
of 1 GeV/nucleon, notwithstanding the low-energy approxi- 
mations. Its low-energy behavior is very similar to that of 
HM, which uses S =0.6 in the whole halo. At energies 
above ~ 50 GeV/nucleon, THMb predicts a significant flat- 
tening for the B/C ratio, which is also suggested by recent 
data from TRACER. Similar conclusions can be drawn for 
THMa. We also expect a similar trend (but less obvious) for 
the p/p ratio, which should be flatter than that estimated in 
Donato & Serpico (201 1). All these features can be tested by 
the upcoming Alpha Magnetic Spectrometer (AMS) 1 data on 
elemental spectra (H to Fe) and nuclear ratios (B/C or p/p). 
We remark that a multi-channel study may be necessary to 
disentangle this effect from other possible processes, see, e.g., 
Tomassetti & Donato (2012). 

The hardening of the B/C ratio is also connected with the 
anisotropy. This can be illustrated for the latitudinal compo- 
nent of the anisotropy, rj z , which is due to the vertical outflow 

1 http : / /www . ams02 .org 



cA(E) h 



(10) 



where z jh = 0.2 (Jones et al. 2001, and references therein). 
Figure 2b shows rj z (E) for the considered models compared 
with direct measurements of the anisotropy amplitude. Due 
to limitations implicit in Eq. 10 as well as possible contribu- 
tions from stochasticity (Blasi & Amato 2012b), the model 
comparison with the data is necessarily qualitative. We fo- 
cus on the comparison among the models, which is more 
instructive. rj z increases with energy as ~ii"o/A (which is 
fixed by the B/C ratio and cannot be arbitrarily changed) so 
that the HM yields r\ z oc E s , whereas THMa and THMb 
predict a natural flattening of r\ z . In particular, THMb pre- 
dicts a much lower anisotropy while keeping a good agree- 
ment with the B/C data at E < 100 GeV/nucleon. Its high- 
energy behavior r\ z oc E 1 ^ 6 is similar to that of Scenario P in 
Vladimirov et al. (201 1), but our model does not fall into any 
of their scenarios. It is interesting noticing that the anisotropy 
may also be reduced at all energies if one accounts for a 
proper radial dependence for K. This is shown in Evoli et al. 
(2012), which assumed a spatial correlation between diffusiv- 
ity and SNR density. In our THMb model, the source spectra 
are fairly soft {v > 2.3) for the DSA requirements, however 
they are not exceedingly soft according to recent calculations 
where the Alfven speed of the scattering centers is accounted 
(Caprioli 2012). On the other hand, THMa gives a depen- 
dence 7] z oc E 1 / 3 at high energies, but requires v w2.2, 
which agrees better with the basic DSA predictions. We re- 
call that state-of-the-art models employ v « 2.4 and S as 1/3 
(Strong et al. 2007) and require large amounts reacceleration 
to match the B/C data. Besides, other B/C+p/p combined 
analyses favor 8 w 1/2 and smaller amounts of reacceleration 
(Di Bernardo et al. 2010). Within our scenario, the observed 
steepness of the low-energy B/C ratio shares the same origin 
of the high-energy hardening of the primary spectra. It can be 
therefore explained, under a purely diffusive picture, why all 
HM-based studies lead to systematically large values for the 
parameter 5. 

4. CONCLUSIONS 

In this Letter, we have proposed a new interpretation for 
the spectral hardening observed in CRs. As shown, it may 
be a consequence of a spatial change of the CR diffusion 
properties in the Galaxy. From this scenario the hardening 
arises naturally as a local effect and vanishes gradually in 
the outer halo, where the CR spectra are also predicted to 
be steeper. This effect must be experienced by all CR nuclei 
as well as by secondary-to-primary ratios. Recent data {e.g., 
Fenni/LALT or TRACER) seem to support this scenario, but 
the predicted spectral upturn is more gradual than that sug- 
gested by PAMELA data. With dedicated analysis of the data 
forthcoming from AMS or Fermi/LAT, our model can be res- 
olutely tested and discriminated against other interpretations. 

Our scenario has remarkable implications for the CR 
physics. With the THMa setup, we have shown that a Kol- 
mogorov diffusion for the inner halo (6 ~ 1/3) is consistent 
with relatively hard source spectra (v ~ 2.2). With the THMb 
setup, we have shown that a very slow diffusion for the inner 
halo (S ~ 1/6) can reconcile the anisotropy with the B/C ra- 
tio. Interestingly, it does not require prohibitively steep source 
spectra as one might expect. Both formulations are based on 
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plain diffusion models where the diffusion coefficient is not 
separable into energy (rigidity) and space terms as usually 
assumed. Their good agreement with data suggests that the 
CR propagation at < 50 Ge V/nucleon might be only moder- 
ately affected by low-energy effects such as reacceleration. 
At the highest energies, we found that the elemental knees 
can be matched using ~ 2.5-10 6 GV of maximum SNR rigid- 
ity, which is attainable by known acceleration mechanisms. 

Further elaborations may be performed using numerical 
methods, in order to introduce a radial dependence of K or 



other effects connected with the Galactic structure. Calcula- 
tions of other observables such as p/p, 10 Be/ 9 Be, or 7-rays 
may lead to a deeper understanding of the effect. We hope 
that our proposal will stimulate further investigations on this 
subject. 

I thank F. Donato for insightful suggestions and B. Bertucci 
for constant encouragement. This work is supported by the 
Italian Space Agency under contract ASI-INFN 1/075/09/0. 
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